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Abstract — The dynamics of many systems nowadays follow not 
only physical laws but also man-made rules. These systems are 
known as discrete event dynamic systems and their performances 
can be accurately evaluated only through simulations. Existing 
studies on simulation-based optimization (SBO) usually assume 
deterministic simulation time for each replication. However, in 
many applications such as evacuation, smoke detection, and 
territory exploration, the simulation time is stochastic due to 
the randomness in the system behavior. We consider the com- 
puting budget allocation for SBO's with stochastic simulation 
time in this paper, which has not been addressed in existing 
literatures to the author's best knowledge. We make the following 
major contribution. The relationship between simulation time 
and performance estimation accuracy is quantified. It is shown 
that when the asymptotic performance is of interest only the 
mean value of individual simulation time matters. Then based 
on the existing optimal computing budget allocation (OCBA) 
method for deterministic simulation time we develop OCBA 
for stochastic simulation time (OCBAS), and show that OCBAS 
is asymptotically optimal. Numerical experiments are used to 
discuss the impact of the variance of simulation time, the impact 
of correlated simulation time and performance estimation, and 
to demonstrate the performance of OCBAS on a smoke detection 
problem in wireless sensor network. The numerical results also 
show that OCBA for deterministic simulation time is robust even 
when the simulation time is stochastic. 

Index Terms — Discrete event dynamic system, simulation- 
based optimization, optimal computing budget allocation. 

I. Introduction 

The dynamics of many systems nowadays follow not only 
physical laws but also man-made rules. These systems are 
known as discrete event dynamic systems (DEDS's). Simu- 
lation is usually the only faithful way to accurately describe 
the dynamics of such systems. The performance optimization 
of these systems then enter the realm of simulation-based 
optimization (SBO). Most existing studies on SBO assume 
deterministic simulation time for each replications. However, 
there exist a large set of DEDS's where the simulation time is 
stochastic. Estimating the evacuation time for a building, the 
smoke detection time of a wireless sensor network, and the 
exploration time of a multi-agent system under a collaborative 
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search policy are such examples. It is of great practical interest 
to allocate the computing budget among designs so that the 
best design can be found with high probability. However, to the 
author's best knowledge, this problem has not been considered 
in existing literatures. 

We consider this important problem in this paper. 
Simulation-based optimization with stochastic simulation time 
is nontrivial due to the following difficulties. First, simulation- 
based performance evaluation. Simulation is usually time- 
consuming, and only provides noisy estimations. In order to 
obtain an accurate performance estimation, one has to run 
simulation by infinite times, which is infeasible in practice. 
Second, discrete variables. Designs usually take discrete and 
finite values. This not only makes traditional gradient-based 
search algorithms not applicable, but also makes the size of the 
search space increase exponentially fast when the system scale 
increases, which is also known as the curse of dimensionality. 
Third, the huge number of computing budget allocations. One 
usually does not have time to explore all the allocations to find 
the optimum. Instead, sequential allocations that can iteratively 
improve their performances are of more practical interest. 
Fourth, stochastic simulation time. Giving the total simulation 
time that is allocated to a design, it is not clear how many 
replications can be finished. Thus the resulting performance 
estimation accuracy is not clear. 

There exist abundant literatures to address the above first 
three difficulties. For example, ranking and selection (R&S) 
procedures are typical procedures for SBO. Bechhofer et al. 
(2, Swisher et al. 0, and Kim and Nelson [3] provided 
excellent review of the R&S works. Chen ||4], Chen et al. 
0, Chen et al. J6), and Chen and Yucesan [7] developed 
the optimal computing budget allocation (OCBA) procedure to 
maximize the probability of correctly selecting the best design 
under a given computing budget. OCBA has been shown to 
perform asymptotically optimally when the simulation time is 
(identically or nonidentically) deterministic. OCBA has been 
extended to tackle the case where the deterministic computing 
time for one simulation replication is different across the al- 
ternatives [8|, to handle multiple objective functions ||9), iflOll . 
simulation-based constraints |11|, opportunity cost [12|, and 
complexity preferences [ 1 3 1 . A comprehensive introduction 
to OCBA is recently available in [14|. Recent good surveys 
on other methods for SBO can be found in lfT5l -H9l. The 
above existing literatures assume deterministic simulation time 
and do not address the aforementioned difficulty of stochastic 
simulation time. 
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In this paper we consider the computing budget allocation 
for SBO with stochastic simulation time and make the fol- 
lowing major contribution. The relationship between the total 
simulation time and the accuracy of performance estimation is 
quantified. It is shown that when the asymptotic performance 
is of interest only the mean value of individual simulation 
time matters. Then based on OCBA for deterministic simu- 
lation time we develop OCBA for stochastic simulation time 
(OCBAS) and show that OCBAS is asymptotically optimal. 
Numerical experiments are used to discuss the impact of the 
variance of simulation time, the impact of correlated simu- 
lation time and performance estimation, and to demonstrate 
the performance of OCBAS on a smoke detection problem in 
wireless sensor network. The numerical results also show that 
OCBA for deterministic simulation time is robust even when 
the simulation time is stochastic. 

The rest of this paper is organized as follows. We present 
the mathematical problem formulation in section HU provide 
the main results in section UTH show the numerical results in 
section IIV1 and briefly conclude in section [V] 

II. Problem Formulation 

Consider a finite set of designs 6 = {1, ...,&}. Let Ji 
be the true performance of design i, which can be accurately 
evaluated only through infinite number of replications 



J, = lim 



E 



where n is the number of replications that are used, 

M£j) = Jt + W i(£]), 

£j represents the randomness in the j-th sample path, and Wi 
has i.i.d. Gaussian distribution N(0, of). Let ti denote the 
time that is consumed by an individual replication of design i. 
We assume the simulation is conducted by a digital computer 
and thus U takes positive integer values and is stochastic. 
Let fi and Fi be the probability density function (PDF) and 
cumulative distribution function (CDF) of ti, respectively. 
Assume that Ji and ti are mutually independent. The case 
when Ji and fcj are correlated will be discussed in section |IV] 
Giving T i7 the number of replications that design i can be 
simulated is stochastic, which is denoted as rii. Then we have 



small to large according to J,;, and denote the best design as 
b. Define the probability of correct selection (PCS) as 



■{Jb 



PCS = Pr < Jb < Ji,i — I, ■■■ ,k and i 



In other words, we are interested in the probability that the 
observed best is the truly best. Now we can mathematically 
formulate the problem as 



Ti,...,T, 



PCS s.t. T * 

i=l 



T. 



where T is the total computing budget. In the following 
discussion we will refer this problem as PI. In other words, 
we are looking for an allocation of the simulation time among 
the designs so that the probability of correct selection is 
maximized. Note that regarding tij as the event occurrence 
time and Ji(£j) as the reward, then {rii(Ti) Ji(Ti)} is a 
renewal reward process [20 . This fact will be used to show 
that Gi can be approximated by a Gaussian (Lemma |3). 



III. Main Results 

In this section, we address problem PI in three steps. 
First, the relationship between the simulation time Ti and the 
distribution of Ji is quantified. Second, an approximation of 
the PCS is provided. Then an approximate version of problem 
PI is formulated and denoted as problem P2. Third, OCBAS is 
developed and is shown to solve P2 asymptotically optimally. 

We start from quantifying the distribution of rn. 

Lemma 1: For any nonnegative integer c, we have 



Pr {rii > c} = 



(F^fr'm, c >i, 

1, c = 0, 



where [/ * g](t) represents the convolution of / and g, i.e., 



[f*g](t) = / f{r)g{t-T)dr- 

J — OO 

r = /*•••*/, 



where tjj is the simulation time for the j-th simulation of 
design i. Assume that an incomplete simulation does not 
output any estimate. When is large, it is reasonable to 
assume that rii > 0. The estimate of Ji is 

1 ni 

■^-E-to 

J=l 

We take the Bayesian viewpoint, which means that the 
estimates Ji , . . . , J^ are given and the true performances 
J\, . . . , Jfc have posterior estimates Jx, . . . , Let gi and Gi 
be the PDF and CDF of Ji, respectively. Sort the designs from 



and f° is the Dirac delta function. 

Proof: When c = 0, we have Prjrii > 0} = 1. When 

c > 1, we have Pr {rii > c} = Pr|X^=i*i,j — ^iV Then 
we prove by induction. 

Step 1. c = 1. We have Pr{n 4 > 1} = Pr{t tl < T,} = 
Fi{T t ) — (Fi * f°)(Ti). 
Step 2. Suppose that we have 



Pr{n 4 > c} = Pr \ <T t \ — (Fi * / c_1 )(T , i ). 

J =1 
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,c+i < Ti 



Step 1 has provided one such example for c = 1. We have 

r c+i 

Pr {n, > c + 1} = Pr \ *<j ^ T * 

= / ''' jE'< ' '< 

f C 1 

x Pr < > ti i — x } dx 

i- J 

= y*Pr{t i , c+1 < 2} - a:} d Prj ^ tij ;< 



E*« 



F i (T i -x)f?(x)dx = (F i *f?)(T i ). 

Combining steps 1 and 2 together, we have 

Pr{n l >c} = (F l *fr 1 )(T l ). 

This completes the proof. ■ 
Lemma 2: 

Fr i^- c >- \ i-F^), c = 0. 

Proof: When c = 0, we have 

Pr {ni - 0} = Pr{t 4 ,i >T t } = l- F t (T,). 
When c > 1, we have 

Pr {rii = c} = Pr {rii > c} — Pr {n^ > c + 1} 

= (F^/r 1 )^)-^*/^) 

- [F i *(/r 1 -/f)](r i ). 

This completes the proof. ■ 
Now we have 

Theorem 1: G*(»)=E^i*(^fe)[^ * (iT* " /f)] (H) 
+G?0) [1 - Fi(Ti)], where $(•) is the CDF of N(0, 1) and 
G?(ie) is the a priori CDF of J;. 
Proof: We have 

Gi(a:) = Pr{ji<z} 

OO 

=J2 Pr < x, rii = c} + G°(x) Pr {m = 0} 

c=l 

OO 

=E Pr {^ ^ a:|ni=jPr{ni=c}+G°(a:)Pr{ni=0}.(l) 

c=l 

Note that 



Pr{j i <a;|n i =j=Pr| £ J,fe)/c < 4=$(j^ 



X J i 



, (2) 



where the first equality follows from the assumption that 
and ti are independent, and the second equality follows from 



the assumption that Ji{^j) are i.i.d. Gaussian. Combine Eqs. 
([T} and (O and Lemma [2] we then have 



This completes the proof. ■ 
Theorem [TJ implies that is not Gaussian. Instead, its CDF 
GAx) is a weighted average of a sequence of Gaussian CDF's 



$ f J 's, which has equal mean values and decreasing 
variances. Note that when T — > oo, Pr{rii = c} is almost zero 
for most values of c except for c s» Ti/fii, where /i; = E[t, 
In this case c is the expected value of m. We have 

Lemma 3: rim<r i _>. 00 — Jij < = $ (x/ (Ti) . 

Proof: Following the elementary renewal theorem 
we have MnciTi^oo E[ni]/7j = 1//^- Thus when T"j — > oo, rij 
also goes to infinity. Then Lemma [3] follows naturally from 
the central limit theorem. ■ 

Lemma [3] implies that when Ti is large, Gi (x) can be rea- 
sonably approximated by N(Ji,afni/Ti). The performance 
of the allocation procedure using this approximation will be 
shown by numerical experiments in section [IV] Following the 
Bonferroni inequality we have 

k 

Pr{j fe < Ji,i= l,...,k and i ^ &}>1-E Pr {^ > ^} ■ 

Following the above analysis Gb(x) and Gi(x) can be approx- 

and $ I , x 7 Ji ) , respectively. 



imated by $ 
Then we have 



x-Jh 



*{*>*} ™f_^-j£°*{-j}*> 

where Sb,i = Jb — J% and of ,- = of^b/T& + afm/Ti. Define 
the approximate probability of correct selection (APCS) as 



APCS = 1 



E 



1 



exp 



Then problem PI can be approximated by 

k 

max APCS s.t. Vt, = T. 

Ti,...,T fc f-f 

4 — 1 

Denote the above problem as problem P2. Replacing rii by 
Ti/ Hi in OCBA J6), we omit the proof and directly present 
the following theorem. 

Theorem 2: Given a total computing time T to be allocated 
to k competing designs whose performances are depicted by 
random variables with means Ji, . . , , J), and finite variances 
a\,...,<j\, and whose individual simulations take random 
time with means fix, . . . and finite variances, as T — > oo, 
the APCS can be asymptotically maximized when 
(1) 



Ti 



I % '* , i,j e {1, . . . , fc} and i ^ j ^ b; (3) 
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Algorithm 1 Optimal computing budget allocation for stochas- 
tic simulation time (OCBAS) 

Step 0: Simulate each design by T time; I <- 0; T[ = = 
■■■It I-,. 

Step 1: If Eti^ > T ' ^op- 
Step 2: Increase the total simulation time by At and 
compute the new budget allocation T[ +1 , 
Theorem [2] 

Step 3: Simulate design i for additional max (0, T^ +1 
time, i = 1, . . . , k; I <— I + 1. Go to step 1. 



T l k +1 using 



T 



(2) 



\ 



k 

E 



(4) 



where 7} is the simulation time allocated to design i, 8b : i = 
Jb — Ji, and Jf, = mini Ji- 

Note that in practice the values of Jj's, erf's, and fii's 
usually are not known a priori, and are replaced by the sample 
means and sample variances, respectively. This gives us the 
sequential computing budget allocation in Algorithm 1, which 
is called OCBA for stochastic simulation time (or OCBAS 
for short). Note that each iteration in OCBA fixes the the total 
number of replications and thus takes stochastic time. But each 
iteration in OCBAS fixes the total simulation time. So the 
number of replications of a design in each iteration becomes 
stochastic. Despite this difference, the total simulation time 
allocated to a design in OCBA and OCBAS are very close. 
This will be demonstrated by the close performances of the 
two methods in the next section. 

IV. Numerical Results 

In this section, we present three groups of numerical ex- 
periments to demonstrate the performance of OCBAS. The 
first group discusses the impact of variance of individual 
simulation time (subsection llV-AI ). The second group discusses 
the impact of the correlation between individual simulation 
time and performance estimation (subsection IIV-B) . The third 
group is a smoke detection problem (subsection HV-Cb . Three 
methods are considered. First, equal allocation (EA), which 
equally allocates the simulation time among the designs. 
Second, OCBA, which iteratively allocates the number of 
replications among the designs [6|. Third, OCBAS, which 
iteratively allocates the simulation time among the designs. 

A. Impact of Variance of Individual Simulation Time 

Consider 10 designs with true performances Ji — i — 1, 
i = 1,...,10. The performance estimation of all the de- 
signs have i.i.d. noise iV(0,6 2 ). The individual simulation 
time of each designs are independent and all have the 
same variance. We consider two types of distributions of 
the simulation time. First, uniform distribution. We conduct 
10 groups of experiments to consider 10 values of vari- 
ances, in which the individual simulation time takes in- 
teger values from [11 — j, 9 + j] with equal probability, 




(b) Simulation time has truncated discrete Gaussian distribution 



Fig. 1. PCS's of EA, OCBA, and OCBAS of subsection HV-Al (averaged 
over 10000 replications). 



j = 1,...,10. Second, truncated discrete Gaussian distri- 
bution. We conduct another 10 groups of experiments, in 
which the individual simulation time of design i satisfies 
Pr{*i = x} <x<f>((x-i + 1.5)/j)-<Z>((x-i + 0.5)/j),x = 
1, . . . , 19; Pr{ij = x} — 0, otherwise. Note that different de- 
signs have different truncated discrete Gaussian distributions. 
Assume that the individual simulation time and performance 
estimation are independent. We apply EA, OCBA (no = 
5, A n = 10, which means that each design is observed by 5 
replications in the beginning and 10 replications are allocated 
among the designs in each iteration afterwards), and OCBAS 
(To = 50, At = 100, which means that each design is 
observed using 50 units of time in the beginning and 100 units 
of simulation time are allocated among the designs in each 
iteration afterwards) under T = 1000, . . . , 10000. Note that in 
each iteration of OCBA we calculate the additional number of 
simulations that are allocated to each design, which may take 
a random simulation time to complete. This is different from 
OCBAS, in which we allocate the simulation time directly. 
The PCS's are estimated using 10000 replications and shown 
in Fig. [T] We make the following remarks. 

Remark 1. When the computing budget increases all three 
methods achieve higher PCS's. This shows that all three 
methods can incrementally improve the PCS's when more 
computing budget is available. This salient feature allows 
sequential allocations of the computing budget, which are 
usually preferred over fixed allocations beforehand in practice. 

Remark 2. OCBAS substantially improves the PCS for a 
given T (or in other words substantially saves the computing 
budget to achieve a given PCS). 
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Fig. 2. PCS's of EA, OCBA, and OCBAS of subsection HV-Bl (averaged 
over 10000 replications). 



Remark 3. OCBA and OCBAS achieve very close per- 
formances. The difference between their PCS's are due to 
randomness of the experiments. This is because the allocation 
procedures of OCBA and OCBAS share the same spirit. 
The difference is that each iteration in OCBA fixes the total 
number of replications and thus takes stochastic time. But each 
iteration in OCBAS fixes the total simulation time. Then the 
total number of replications of a design becomes stochastic. 

Remark 4. When the variance of individual simulation 
time increases, the PCS's of OCBA and OCBAS do not 
change much. This is because each design is observed by 
more and more times when T increases. So the variance of 
the simulation time does not significantly contribute to the 
performance estimation when the total simulation time is large, 
as shown in Lemma [3] Note that when the total simulation 
time is small, the variance of the simulation time matters. For 
example, when T < 10 (the mean value of each individual 
simulation), a larger variance allows to complete a simulation 
within T with a larger probability. When the T > 10, this 
impact of the variance reduces fast. 



B. Correlated Simulation Time and Performance Estimation 

Consider 10 designs with true performances Jj = i— 1, i = 

1. . . . , 10. The performance estimation of all the designs have 

1.1. d. noise N(Q, 6 2 ). The individual simulation time U of 
design i takes values of 5 and 15 with equal probability and 
is correlated with its performance observation noise Wi as 
follows. If Wi > 0, then t.- L — 15 with probability (w.p.) p 
and ti — 5 w.p. 1 — p, where < p < 1 is a given constant. If 
Wi < 0, then ti — 15 w.p. 1 — p and ti — 5 w.p. p. The value 
of p indicates the correlation between ti and Wi. In particular, 
p = means that fj and Wi are purely negatively correlated; 
p = 1 means that ti and Wi are purely positively correlated; 
and p = 0.5 means that ti and Wi are independent. For 
p = 0, 0.1, . . . , 1.0 and T = 1000, . . . , 10000, we apply EA, 
OCBA (n = 5, A„ = 10), and OCBAS (T = 50, A T = 100) 
and estimate the PCS's by 10000 replications (shown in Fig. 
01. Remarks 1-3 also hold in this case. We can also see that 
the correlation p does not affect PCS much. 




Fig. 3. A smoke detection problem in wireless sensor network. 



C. A Smoke Detection Problem 

We compare three methods on a smoke detection problem 
in wireless sensor network. Consider an area of interest (Aol) 
with unit size as shown in Fig. [3] which is discretized into 
11 x 11 grids. A fire may be set at any point on the grid 
inside the Aol with equal probability. Once a fire is set on, 
it generates a smoke particle within each time slot. A smoke 
particle may walk to a neighboring grid in each time slot. 
There are at most four such neighboring grids corresponding 
to four directions. The probability to walk to one of the four 
grids is proportional to its distance to the fire source, i.e., 

Pr{x t+1 = x t + l,y t+1 = y t } cx d((x t + 1, y t ), (x , y )), 

Pr{x t+ x=x t -l,yt+i = yt} oc d((x t - 1, y t ), (x , yo)), 

Pr{x t+ i = x t ,y t+1 = y t + 1} cx d((x t , y t + 1), (x , y )), 

Pr{x 4+ i = x t ,y t +i = yt - 1} oc d((x t , y t - 1), (x Q , y Q )), 

where (xo,yo) represents the position of the fire source and 
d(-, ■) represents the distance between two positions. Once a 
particle walks to the boundary of Aol, it is bounced back. 
There are 3 sensors that can be allocated to the 9 positions 
marked by circles in Fig. [3] Once a smoke particle arrives 
at any of the three sensors, it is detected. The question is 
how to allocate the sensors to minimize the average detection 
time. It is easy to show that there are 84 allocations in total. 
Considering the symmetries, only 16 allocations need to be 
considered. The response time of the 16 designs are evaluated 
by 100000 independent replications and shown in Table [I] 
where the designs are represented by the positions of the three 
sensors. Note that in each simulation the response time takes 
integer values. But the mean values of the response time take 
positive real numbers. As an example, we show the probability 
mass function of the response time of the first design (design 
1,2,3) in Fig. |4j which is estimated by 100000 independent 
replications. Note that in this example we have Ji = ti, i.e., 
the performance estimation and individual simulation time are 
the same. This violates the assumption used in Theorem [2] 

The probability of correct selection of the three methods 
for T = 1 x 10 4 , 2 x 10 4 , . . . , 1 x 10 5 are evaluated using 
10000 independent replications and shown in Fig. [5] Since 
each individual simulation of a design takes about 10 units of 
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Fig. 4. The probability mass function of the response time of the first design, 
estimated by 100000 replications. 



TABLE I 

True performances of the designs 



Index 


Design 


Response time 


Index 


Design 


Response time 


1 


1,2,3 


11.5989 


9 


1,3,7 


8.6777 


2 


1,2,4 


10.6383 


10 


1,3,8 


7.6482 


3 


1,2,5 


9.3776 


11 


1,5,6 


8.0903 


4 


1,2,6 


9.3353 


12 


1,5,9 


8.1355 


5 


1,2,7 


9.7781 


13 


1,6,8 


7.4699 


6 


1,2,8 


8.2390 


14 


2,4,5 


8.5127 


7 


1,2,9 


8.7794 


15 


2,4,6 


7.6968 


8 


1,3,5 


8.6344 


16 


2,5,8 


7.7671 



time (as shown in Table IJi to make a fair comparison between 
OCBA and OCBAS, we use the following parameter settings. 
In OCBA, let n = 20 and A„ = 10. In OCBAS, let T = 
200 and Ay = 100. Remarks 1-3 also hold in this case. We 
can see that OCBAS works well even when the performance 
estimation and individual simulation time are correlated. 



V. Conclusion 

In this paper, we consider the computing budget allocation 
for SBOs with stochastic simulation time and develop OCBAS 
to provide the allocation asymptotically optimally. The per- 
formance of OCBAS is demonstrated through two groups of 
academic examples and a smoke detection problem in wireless 
sensor network. The numerical results also show that OCBA 
for deterministic simulation time is robust even when the sim- 
ulation time is stochastic. Note that the asymptotically optimal 
allocation of computing budget only depends on the mean 
value of the stochastic simulation time. Note that though we 
assume the performance estimate Ji and the simulation time 
for an individual replication ti are independent in section [HI] 
the numerical results in section[IV]show that OCBAS performs 
well when Ji and i, are correlated. Note that Lemma [3] shows 
that the performance estimator under stochastic simulation 
time can be well approximated by the performance estimator 
under deterministic simulation time. Replacing rii by Ti/ Hi, 
OCBAS can be obtained straightforwardly from OCBA. Using 
Lemma [3] it is possible to extend OCBAS to handle multiple 
objective functions, simulation-based constraints, opportunity 
cost, and complexity preferences, following its according 
extensions in OCBA. That will be important future work. 
Note that when N parallel computers are available, the total 
computing budget we can use will be improved from T to NT. 
Both OCBA and OCBAS can be extended to this situation. 
But if each computer can simulate only a specific design, i.e., 
N = k, we usually have a constraint on the decision making 
time that is max^ Ti. How to allocate computing budget 
according to this constraint is an interesting future research 
topic. We hope this work brings insights on addressing SBOs 
with stochastic simulation time in general. 
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Fig. 5. PCS's of EA, OCBA, and OCBAS of subsection ITV-Cl (averaged 
over 10000 replications). 
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